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Nonlinear Patterns in Urban Crime: Hotspots, Bifurcations, and Suppression* 
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Abstract. We present a weakly nonlinear analysis of our recently developed model for the formation of crime 
patterns. Using a perturbative approach, we find amplitude equations that govern the development 
of crime “hotspot” patterns in our system in both the one-dimensional (ID) and two-dimensional 
(2D) cases. In addition to the supercritical spots already shown to exist in our previous work, 
we prove here the existence of subcritical hotspots that arise via subcritical pitchfork bifurcations 
or transcritical bifurcations, depending on the geometry. We present numerical results that both 
validate our analytical findings and conhrm the existence of these subcritical hotspots as stable 
states. Finally, we examine the differences between these two types of hotspots with regard to 
attempted hotspot suppression, referencing the varying levels of success such attempts have had in 
real world scenarios. 
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1. Introduction. The study of pattern formation in physical and mathematical systems 
has a long and interesting history. This general subject area is also quite diverse, examining 
biological (see, as a small sample, [37, 28, 27]), geological [33, 10, 2], and even sociological 
systems [30, 14], to name but a few. Though these various subjects and systems may seem 
completely unrelated, the mathematics describing the patterns in each are surprisingly similar. 
Consequently, a robust, powerful, and universal set of mathematical tools has been developed 
to study such systems, and the employment of these tools can lead to better understanding 
of pattern forming systems, regardless of their specific nature. 

Recently, we set forth to develop a mathematical model to describe the spatio-temporal 
patterns of urban crime [35]. Using well-known criminological ideas regarding the way in which 
criminal events affect future crime risk in a location, and the way in which risk can spread from 
one area to another [18, 19, 20, 1], we constructed a model consisting of two coupled, nonlinear 
PDEs that may describe the formation and dynamics of crime “hotspots”—spatio-temporal 
clusters of high crime. Using a simple linear stability analysis of our model, we found that the 
homogeneous system can be unstable to disturbances of specific wavenumbers under certain 
parameter regimes, leading to hotspot formation. However, our previous work stopped there, 
with no investigation of the possibility of hotspots outside of this linearly unstable regime. This 
paper addresses this possibility by performing a weakly nonlinear analysis on our system and 
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developing amplitude equations for the model; this is a detailed follow-up to our paper [34], 
which presents only a few of the qualitative results of such an analysis. By investigating 
the possible bifurcations in the steady state solutions of our system both analytically and 
numerically, we indeed find that stable, “large” amplitude hotspots may exist even in the 
linearly stable regime. 

The fact that these subcritical hotspots exist within our system is especially interest¬ 
ing when attempting to understand the outcome of hotspot suppression, typically by police 
executing a strategy known as hotspot policing, which has become dominant over the past 
two decades [6, 4, 5, 41, 40]. Recognizing that crime tends to form dense clusters in space 
and time, leaving some areas with little or no crime problem, police routinely target their 
limited resources at those locations experiencing high crime. That hotspot policing would 
be an improvement over random patrol is uncontroversial; it has been well known since the 
1970s that random patrol has little measurable effect on crime [23]. However, questions have 
been raised about whether hotspot policing leads to lasting hotspot reductions, or simply the 
displacement of hotspots from one area to another [29, 3, 7]. The present research provides 
a formal theoretical foundation for understanding different potential outcomes from hotspot 
policing in relation to the classification of hotspots as either supercritical or subcritical [34] . 

The remainder of the paper is organized as follows. In section 2, we give a brief introduction 
to our crime model and the major results found in [35]. In section 3, we perform a weakly 
nonlinear analysis of our system in both the one-dimensional (ID) and two-dimensional (2D) 
cases, deriving some analytical results for the amplitude equations and bifurcations governing 
the hotspots exhibited by the system. In section 4, we compare these analytical results to 
numerical solutions. Finally, in section 5, we explore the possible results of hotspot suppression 
qualitatively and numerically using both the continuum and discrete models. 

2. Background. We begin by reviewing the results of [35]. First, we developed an agent- 
based model of criminal activity that aims to reproduce the known phenomena of repeat and 
near-repeat victimization [18, 19, 20, 1], whereby crime risk becomes elevated in an area and 
its surroundings following an initial event there. This model couples the dynamics of moving, 
offending criminals on a 2D lattice (with lattice spacing i) with an underlying scalar field 
H(x, t) that we refer to as the attractiveness. As the name implies, the attractiveness field is a 
measure of how desirable any given location on the lattice x = (i, j) is as a target for criminal 
activity, with the numerical value of the field giving the stochastic rate of offending for the 
n(x, t) criminals at that location. The model evolves in discrete time, using a timestep 6t, 
and during each timestep criminals may victimize their current location with probability 

(2.1) p,(x,t) = l-e-^("’‘)'5'. 

If the criminal does in fact choose to commit a crime during this timestep, he is removed 
from the lattice. If, on the other hand, he does not commit a crime during this timestep, he 
will instead move to one of the four lattice points adjacent to his current location, selecting a 
particular neighbor x' with probability 

(2.2) p„(x'.t;x)= , 

x"~x 
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Figure 1. Example output from the discrete system. These colormaps display high A in red, A in green, 
and low A in blue to purple. In (a) is an example of stationary hotspots, where once hotspots form at a point, 
they tend to stay there indefinitely; in (b) is an example of transitory hotspots, where hotspots do not typically 
last forever, and will move about, deform, or disappear over time; and in (c) is an example of no hotspots, 
where no large spots are ever observed. 


where the notation x" ~ x indicates all of the sites neighboring site x. In this way, criminals 
actively seek out areas of high A, where they are more likely to commit crimes. 

The attractiveness held is composed of a static component AP (referred to as the “baseline” 
attractiveness) and a dynamic component i?(x, t), such that ^(x, t) = i?(x, t)+A^. After all of 
the criminal activity for the round is completed, the dynamic component of the attractiveness 
held is updated through the following mechanisms, meant to model repeat and near-repeat 
victimization. First, the dynamic attractiveness spreads spatially via a weighted averaging 
procedure between each site and its four neighbors. Next, it exhibits an exponential decay 
with rate oo. Finally, for each criminal event that occurred at x, the attractiveness is increased 
there by an amount 9. So, if the number of events at site x during the current timestep is 
given by E{x), then the attractiveness at the beginning of the next timestep will be 


(2.3) 


A(x, t + St) 


(1 - r/)B(x, t) + ^Yl 

x'~x 


{l-ujSt) + dE{^)+A°. 


Note that, by dehnition, i] < 1. As a hnal step of the simulation, new criminals are generated 
on each point of the lattice at a rate F. 

The discrete model thus described contains a number of parameters, and depending upon 
the choice of these parameters, the system may exhibit three general types of behavior: sta¬ 
tionary (fixed in space) crime hotspots, transitory (moving about in space or appearing and 
disappearing in time) hotspots, or no hotspots at all; these three cases are illustrated in 
Figure 1. 

Second, we derived a continuum limit of the discrete model by converting the criminals 
into a number density p{x, t), taking expectation values for all probabilistic events, and letting 
£, —>• 0 with the constraint £‘^/6t = D, a diffusion coefficient ([35] presents all of the algebraic 

details). This hydrodynamic limit results in our PDF model, which is the major focus of this 
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paper and can be written in the dimensionless form 

, /tO , „ 


(2.4) 

dA 2 4 

a 

(2.5) 

dp ^ 


A 


- pA + A-A°, 


where rj is the same as in the discrete model, and the two remaining parameters A and A^ can 
be found from various combinations of the six remaining parameters present in the discrete 
model. Inspecting these equations, we see that crimes occur locally at rate pA, and each such 
crime causes A to increase. In addition, A diffuses with dimensionless diffusion coefficient p 
(< 1) and decays exponentially to the baseline value A^. Criminals exhibit diffusive motion 
with an advective bias up gradients of In A. Finally, criminals are subtracted from the system 
when they commit a crime and are added back at a constant rate A — A^} These equations 
exhibit a general reaction-diffusion form and are similar to models of other biological systems, 
such as the Keller-Segel chemotaxis model, which are well studied in the literature (see, for 
example, [22, 16, 39, 9, 13, 24, 36, 15, 8, 31]). 

The continuum system described by (2.4) and (2.5) may display two of the three behaviors 
from the discrete system: stationary hotspots or no hotspots. We believe that transitory 
hotspots are not seen in this continuum approximation since they are the result of statistical 
noise that is removed by considering only expectation values in our limit. We showed that 
the formation of hotspots in this system may arise as a result of a linear instability of the 
homogeneous steady state 

— A^ 

(2.6) A{x,t) = A, p(x,t) = p=l-^ 

toward perturbations of certain wavenumbers k, and that the dispersion relation could be 
written as 


(2.7) u'(k) — — [1 + j 4 — p + |k|^(l + ?])] /2 

+ \J[l + A-'p+ |k|2(l r/)]^ /4 - (r/jkj^ - (3p - pA - l)|k|2 + A). 

The instability criterion, therefore, could be written as 

( 2 . 8 ) ^0 < ^0 ^ 2 - _ 1^-2 _ 

In other words, if the baseline attractiveness is less than some critical value A^, then the 
homogeneous state will be linearly unstable (exhibit some modes with a positive a). Finally, 
we showed that the maximally unstable mode k^ax is given by 

(2.9) jkmaxl^ = (1 - ^)/(1 - ?/) - P(5 - p)/{l - pf 

+ ?/(l + ■nf'P [(^(3 -p)-2){l-p) + 2p(3 - p)] /p{l - pf. 

^The choice of the notation A here is due to the fact that, at steady state, this quantity is indeed the 
spatially averaged value of j4(x), regardless of the other parameters or whether hotspots are displayed or not. 
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Note for future reference that, when the maximally growing mode can be greatly 

simplified to 


( 2 . 10 ) = | k *|2 

3. Weakly nonlinear analysis. Our goal now is to more deeply examine the continuum 
system of (2.4) and (2.5) and to move beyond the simple linear stability analysis outlined 
above, thus providing the technical details of the qualitative results presented in [34], and for 
an even wider range of possibilities. We will accomplish this by means of a weakly nonlinear 
analysis, using a standard perturbative expansion approach to derive amplitude equations for 
our system [12, 38, 17, 25, 11]. 

We begin by considering a parameter regime such that the homogeneous state is linearly 
unstable (or stable). Choosing as our control parameter (as suggested by our form of the 
stability criterion given in (2.8)), we define a new parameter e via the equation 

(3.1) 41° = 71° - eA, 



such that the homogeneous state will be linearly unstable for positive e and linearly stable 
for negative e, as indicated by the linear stability criterion in (2.8) above. For the remainder 
of our analysis we will assume that e is small, but we note here that in theory e can take on 
any value between some emin (where ^4^ = A, the most it could ever physically be) and Cmax 
(where ^4^ = 0, the least it could ever physically be). These two values are given by 


(3.2) 

(3.3) 



with the difference between these always being 1. We further point out that Cmax becomes 
negative for any rjA = r]k1 > ^/3 — 1, meaning that above this threshold the homogeneous 
state is incapable of being linearly unstable and our analysis is invalid in this regime. 

Returning again to the results of the linear stability analysis, when we substitute (3.1) 
into (2.7) and expand for small e, we find that the growth rate for the k* mode is given by 


(3.4) 


iT(k*) = (T*e + O(e^), 


where 

(1 + ?7lk*12) [2r/+ r/lk*12(3 - r?)] ■ 

Using this result, we see that we can define a new, slow time variable T = \e\t that describes 
the dynamics of the system when near the stability transition; this means that the dt in (2.4) 
and (2.5) becomes Jej^r. We use jej here to make our future results valid regardless of the 
sign of e, though this means that we must rewrite (3.1) as 

(3.6) - sign(e)lelA. 
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At this point we define a new spatial variable x= |k^,|x (but continue to refer to x as x in the 
future for notational simplicity) and rewrite (2.4) and (2.5) as 


(3.7) 


BA 

|e|^ = - A + - sign(e)r7|k*|'^|€| + pA, 


dp 


(3.8) |e|^ = |k,|2v 


Vp - ^VA 
A 


- pA + p\k^\‘^ - A° + sign(e 
Next, we express A and p as expansions in our small parameter of the form 


k*|‘^|e| 


(3.9) 

(3.10) 


A(x,r) = A + ^|e|“^'A(j')(x,r), 

i=i 

aO °° 

p{^,T) = 1 - ^ + ^|e|"V^^^(x,r), 


where a is a rational number that will depend upon the specific geometry in which we are 
interested; the reasoning behind the choice of a will be presented with each geometry we 
consider. We substitute these expansions into our differential equations and then separate the 
resulting equations by powers of |e|. We note that, upon doing this, (3.7) can be used to simply 
solve for a given p^^^ (x, T) algebraically in terms of lower order p^^ ^ (x, T) and A^^ ) (x, T) and 
their derivatives, and that this result can then be substituted into (3.8). This leaves a series of 
fourth order differential equations to be solved that involve only the various A^^\'x,T), each 
of which is of the form 


(3.11) 


(v^ + i)2A(^)(x,r) = fj [A(^)(x,r) 


where fj is a possibly nonlinear function. Regardless of a, the first of these equations is always 
(3.12) (V2 + l)2A(^)(x,r) = 0. 


3.1. ID case. In this geometry, we restrict our solution to a domain x G [0, L], where 
L = 2mT for some integer n > 0, and impose periodic boundary conditions for both A(x, T) 
and p{x,T). The solution to (3.12) in this geometry and for these boundary conditions is 

(3.13) A(i)(x,r) = P(r)e“+ C.C., 


where P{T) is the amplitude, which at this point is simply an integration constant, and “c.c.” 
denotes the complex conjugate. Due to the inversion symmetry of this solution (P —>• —P is 
physically the same but just shifted), we expect a pitchfork bifurcation to occur here; we use 
a = 1/2 to reflect this. The first interesting equation therefore occurs at order |e|: 

(3.14) (V^ + 1)^ A(2)(x,r) = [P(r)2e2“ +C.C.] . 

The particular solution to this equation, which is all we are after, is 

(3.15) T) = [P(r)2e2“ + c.c.] . 


Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 









468 


M. B. SHORT, A. L. BERTOZZI, AND P. J. BRANTINGHAM 


At order we find the equation 

(3.16) (V2 + 1)' T) = / 3 ,i [P{T)-rj, K] e“ + / 3,3 [P{T)-,rj, h] + c.c.; 

we do not reproduce the full expressions for /sj here for the sake of simplicity. Note that (3.16) 
contains a secular term oc e“. In order for the particular solution of (3.16) to fit the periodic 
boundary conditions chosen, this secular term must vanish, meaning that / 3 q [P{T);r], fc*] = 0. 
Upon enforcing this constraint, rescaling T back to t, and letting \e\^^‘^P{t) = Q{t), we find 
the amplitude equation 

(3.17) Q = a^eQ - Ci{'q,k^)\Q\^Q, 


where 

(3.18) 


Ci{r],k^) 


—8 + 56rjk‘^ — 31r]‘^kf — Srj^k^ 
3ri^kl [2r? + r?A:2(3-r?)] 


and (7* is given by (3.5) above. As expected, this is the standard form for a dynamical 
system exhibiting a pitchfork bifurcation, with the distinction between a supercritical and 
subcritical bifurcation determined by the sign of Ci. Upon inspection, it is found that Ci will 
be negative for any rjkl < 0.157 (indicating a subcritical pitchfork bifurcation) and positive 
otherwise (indicating a supercritical pitchfork bifurcation). The steady state value Qs is either 
zero (the homogeneous case) or given by 


(3.19) 


Qs = P 


cr*e 


Ci{r},k^y 


Finally, our solution for Q is valid only to order |e|, so our solution for A{x,T) is also valid 
only to this order and is given by 


(3.20) 


A{x, t) = A + Q{t)P'^ 


4 (1 - r^kt) 




+ c.c. 


One can in general continue the expansion up to higher orders in e by defining subsequent 
slow timescales Tj for j > 2, each of which will modify dt by adding a term \e\^dTj- One 
then continues with the above results and eliminates the secular terms at higher orders in the 
expansion, with the net result being amplitude equations that govern the various Pxy For 
example, the next order amplitude equation for the ID case takes the form 

(3.21) Q = (T*e [1 + eai(r 7 . A:*)] Q - Ci{r], /c*) [1 + ea 2 {r], /c*)] \QfQ - C 2 {r], k^lQ^Q, 


where oi, 02 , and C 2 are the new corrections that arise as we move to the higher order (the 
exact formulas for these expressions are unimportant here). We will refer to this higher order 
amplitude equation later when discussing numerical simulations to help explain some of the 
results seen there. 
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3.2. 2D, radially symmetric case. We now consider solutions on a disk r S [0,i?] with 
R = where is the nth root of the Bessel function Ji(r’); we enforce Neumann 

conditions on the boundary edge. For these boundary conditions in this geometry, the solution 
to (3.12) is 

(3.22) A(i)(r,r) =P(r)Jo(r). 

In this regime, therefore, the inversion symmetry of the ID case is broken {P —)• —P is 
physically different here), so we expect a transcritical bifurcation rather than a pitchfork 
bifurcation. Hence, we choose a = 1, and the first interesting equation in our system is 
proportional to |ep: 


(3.23) 


(y'^ + if T) 

_ 9r]ky{T) sign(e) - (l + r/A;^) [2r] + r]k‘^{3 - rj)] Pt(T) 


3iffki 


Mr) 


+ 


[4(r) - Jf(r)] . 


As before, we will need to eliminate any secular term proportional to Jo{r) on the right- 
hand side of (3.23) so that our solution will respect the Neumann boundary conditions im¬ 
posed. In order to do so, we take advantage of the fact that the Bessel functions can be used 
as an orthogonal basis for expanding other functions, so we are free to write the Jq (r) — (r) 

portion on the right as a sum of Bessel functions to the first power, one of which will be Joir)- 
With the definition that 

2 f rJo(r) [J^{r) - jf{r)] dr 

Jo 


(3.24) 


Q = 


R^MR) 


we see that setting the secular term to zero (and rescaling T to t and letting |e|P(t) = Q{t)) 
is equivalent to the amplitude equation 


(3.25) 


Q = o-^eQ + C3{r], k^)Q‘^ 


where 

(3.26) 


(73(77, fc*) 


6 g(l - r]kl) 
k‘^[2'n + - rj)]' 


As expected, we find that in the 2D, radially symmetric case, our system will undergo a 
transcritical bifurcation. Interestingly, the constant 6*3 will always be positive for any value 
of rjk)^ for which the above analysis is valid. That is, 6*3 would only be negative if rjk"^ > 1, 
but the maximum value of rjk'l for which linear instability is at all possible (which must be 
the case for our analysis to work) is rjkl = \/3 — 1 < 1, as discussed previously. Hence, there is 
really only one qualitatively distinct bifurcation diagram in this case. The steady state value 
Qs in this case is either zero (the homogeneous case) or 


(3.27) 


Qs = 


cr*e 

C3{r],h)' 
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As in the ID case, our amplitude equation is valid only to order |e|, so our equation for A(r, t) 
in this case is 


(3.28) A{r,t) = A + Q{t)Jo{r). 

As stated above, in this geometry there is a physical difference between positive Q and 
negative Q solutions, with the former corresponding to a solution that exhibits a bump in A 
at the origin (hereafter referred to as the “bump solution”) and the latter corresponding to 
a solution that has a ring of high A at the outer edge of the domain (hereafter referred to 
as the “ring solution”). Our theoretical results from (3.27) state that the steady state bump 
solution will exist only for negative e and that it will be unstable, and that the steady state 
ring solution will exist only for positive e and be stable. 

3.3. Fully 2D case. For the fully 2D system, the first order equation (3.12) (along with 
suitable boundary conditions) admits solutions of the form 

N 

(3.29) A(i)(x,r) = + C.C. 

i=i 

for any N, so long as |qj| = 1. We will limit our discussion here to the simple cases of rolls, 
squares, and hexagons, however, as they display simple periodicity. Rolls {N = 1), though, 
are just the ID patterns discussed previously extended into a second dimension, so we need 
not perform any further analysis on them here. 

3.3.1. Squares. We begin our analysis with squares: = 2, qi = x, q 2 = y, a domain 
X G [0,Lx], y G [0,Ly], where Lx = 2mT and Ly = 2rmr for some integers n,m > 0, and 
periodic boundary conditions for both A{x, T) and p{x, T). The first order solution is therefore 

(3.30) A(^)(x,r) = Pi(T)e“ + P2(R)e*^ + c.c. 

Equation (3.30) displays the same inversion symmetry as the ID case; hence we expect a 
pitchfork bifurcation here as well. Consequently, the mathematics in this case follow almost 
exactly as in the ID case above, so we omit the small details. In the end, we arrive at the 
second order solution 


(3.31) A(2)(x,r) = 


4 (1 - 


Pi{Tfe^^^ + P2{Tfe^^y 

+ 9Pi{T)P2{T)p‘^^+y'> + 9Pi(r)P2(R)*e*("^"^) + c.c. 


and amplitude equation for Qi (as in ID, Qi = y^Pi) 

(3.32) Qi = cr*eQi - C'i(?7, fc*)|Qi|^Qi + Cs{ri,h)\Q2\‘^Qi- 

Here, Ci{rj,k^) is as given in (3.18) above and 


(3.33) 


6 (8 - 8r]k‘^ - 9r]‘^kf + Srj^k^) 

y‘^kl[2rj + rjkU3-v)] ’ 
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the amplitude equation for is symmetric with (3.32), with the subscripts 1 and 2 switched. 

The steady states for squares (Qi = Q 2 = Qs) are therefore Qs = 0 or 

^ Ci{rj,h)-Cs{ri,h)- 

For all valid values of r/k'^, Cs > 0, and for all values of < 0.7, Cs > ICil. Hence, squares 
almost always developed through a subcritical pitchfork bifurcation from the homogeneous 
state, only developing through a supercritical pitchfork over the narrow range 0.7 < i]k‘l < 
^/3 — 1. In either case, the roll steady state is unstable to squares. This is because the roll 
steady state is simply the square system with Qi = Qg from (3.19) and <52 = 0. This fixed 
point is always unstable along the direction of increasing < 52 , however, since the eigenvalue in 
this direction is A 2 = o'*e(l + Cs/Ci), which is always positive. That is, either Ci,e > 0 (a 
supercritical roll), where A 2 is clearly positive since Cs > 0, or Ci, e < 0 (a subcritical roll), 
where A 2 is positive since Cs > |Ci|. 

3.3.2. Hexagons. The final patterns we will examine are hexagons: = 3, qi = x, 

q 2 = —+ ^y, qs = —— ^y, a domain x G [0,Lx], y G [0,Ly], where = dnvr and 
Ly = Armrjy/S for some integers n, m > 0, and periodic boundary conditions for both A{x, T) 
and p{x,T). The first order solution is therefore 

(3.35) H(^)(x,r) = Pi(r)e“ + P2(T)e*(-^/2+V3j//2) ^ p^i^j.y{-x/2-V3y/2) p ^ ^ 


The lack of inversion symmetry in this case, along with the resonance qi + q 2 + qs = 0, will 
generally cause our amplitude equations to display a quadratic nonlinearity, as seen in the 
radially symmetric system above, leading to a transcritical bifurcation. Correspondingly, we 
choose a = 1 in this case and follow a derivation similar to the radially symmetric case. We 
find that the amplitude equation for Qi (recall that Qi = |e|Pi in this case) is given by 


(3.36) 
where 

(3.37) 


Qi = a^eQi + Cniv, h)Q2Ql, 


CH{r],h) 


6(1 - r?/cg) 
fcj [2r] + r]k^{3 - r/)] ’ 


the amplitude equations for Q 2 and Qs are given by cyclic permutations of the indices 1, 2, 
and 3 in (3.36). 

The valid hexagonal steady states {Qi = Q 2 = Qs = Qs) are Qs = 0 or 


(3.38) 


Qs = 


(T*e 

CHir],k*)' 


Note that Ch is positive for all valid values of r/fej and is in fact just Cs (see (3.26)) without 
the factor of q, meaning that the bifurcation diagram for hexagons will be quite qualitatively 
similar to that of the radially symmetric state for very small |e| and amplitudes. Therefore, 
and for future reference, we will continue to denote the positive Q solutions as “bumps” and 
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(a) (b) 

Figure 2. ID system with rjk^ — 0.4. In (a) is a bifurcation diagram for the system where dashed lines 
represent unstable branches and solid lines are stable branches; numerical results are in black, and analytic 
results are in red (with circles). We find very good agreement for essentially all t values. In (b) are plots of 
the numeric (black) and analytic (red with circles) solutions for Aamp{t) with e = 0.01 and Q{0) = 0.01; there 
is very good agreement here as well. 

the negative Q solutions as “rings” even in the hexagonal geometry. However, unlike the 
radially symmetric system, the hexagonal ring solution is not stable but is instead a saddle 
point, meaning that the only stable hexagonal steady state predicted by this approximation 
is Qs = 0 for e < 0. However, the quadratic nonlinearity of the hexagons dominates over 
the cubic nonlinearity of the rolls and squares near onset; hence, hexagons are the preferred 
pattern for our model at small |e|. This is not surprising, as it is a generic feature of systems 
that lack inversion symmetry [12], as ours does. 

4. Numerical results. As a verification of our analytical results above, we numerically 
solve our model system in various geometries. For the dynamical system, we use a fully 
implicit Newton-Raphson-based solver; for the steady state solutions, we use a Newton- 
Raphson-based relaxation method. For each case, we look at a quantity we will refer to as 
simply the “amplitude” of A, which we define as 

(4.1) Aamp(t) = ^ 

where D is the domain of the simulation, and |D| is its size. Our measure is, therefore, 
essentially a root mean square (RMS) measure of the attractiveness field’s deviation from the 
homogeneous steady state. Finally, for all simulations in this section, we employ Neumann 
boundary conditions. 

4.1. ID case. In this geometry, our domain B is x € [0,7r/A:*]. The first case we explore 
is a supercritical system, in which rjk'^ = 0.4 {rj = 0.1 and fe* = 2). The two plots in Figure 2 
summarize the results here. Figure 2(a) shows a bifurcation diagram for our system as derived 
by computing the steady state value Aamp(oo) as a function of e, plotting both the analytical 
and numerical results. We find that there is good agreement in this case for essentially all 
e values. Figure 2(b) plots the analytic and numerical solutions for Aamp(0 using e = 0.01 
and (5(0) = 0.01; there is very good agreement here as well. 
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(a) 


(b) 


Figure 3. ID system with rjk^ =0.1. In (a) is a bifurcation diagram for the system where dashed lines 
represent unstable branches and solid lines represent stable branches; numerical results are in black, analytic 
results from (3.19) are in red (with circles), and a higher order analytic solution as in (3.21) is in blue (with 
squares). There is good agreement between the numerics and both analytic solutions for smaller e values along 
the unstable branch, but only the higher order analytic solution predicts the existence of the large amplitude 
branch seen in the numerics. In (b) are plots of the numeric (black) and analytic (red circles and blue squares) 
solutions for Aamp(t) with e = — 0.001 and varying < 3 ( 0 ); the dashed line represents the analytic steady state 
value for this e. The lower line corresponds to Q(0) = 0.0028, with the analytic solution from (3.19) in red 
(with circles); there is good agreement here between the two. The upper line corresponds to Q{0) — 0.0032, 
with the higher order analytic solution in blue (with squares). The agreement between these two is reasonable, 
though the analytic solution predicts a higher steady state value than the numerics. 


The next case we explore is a subcritical system, in which r]k‘^ = 0.1 {r] = 0.1 and A;* = 1), 
with the results shown in Figure 3. Referring to Figure 3(a), the numerical solutions (black) 
display the small amplitude, unstable branch predicted by the theory above, and the numerics 
match the theory (red with circles) well at small e. However, there is also a stable, large 
amplitude branch in the bifurcation diagram that is not predicted by the theory above. As 
alluded to before, however, if we continue our analytic solution to the next higher order in e as 
in (3.21) (blue with squares), we can predict the location of the secondary bifurcation where 
the upper and lower branches meet. Note, however, that this higher order amplitude equation 
is not necessarily valid along the upper branch seen in the numerics, which explains the 
substantial deviation seen there. Figure 3(b) shows the evolution of A ^m p(t) using e = —0.001 
and two different values for <5(0). The first value is <5(0) = 0.0028, which is just slightly below 
the unstable branch, so we expect our analytic results above (red with circles) to be close to 
the numerical results (black). However, the second Q(0) is 0.0032, which is slightly above the 
unstable branch, so our results above cannot be used. Instead, we compare with the higher 
order analytic result (blue with squares), and find reasonably good agreement until t « 800. 
After this time the validity of (3.21) is clearly lost, as the true solution exhibits oscillatory 
behavior before settling down to the steady state, which is something our amplitude equation 
could not predict. 

4.2. 2D, radially symmetric case. In this geometry, D is r G [0,/3iq/A:*]. Furthermore, 
we have chosen a convention whereby bump solutions are shown with a positive Aamp value, 
while ring solutions are shown with a negative one. Simulations use rjk'^ = 0.2 (77 = 0.01 and 
/c* = 2\/5), and the results are shown in Figure 4. In Figure 4(a), we see our bifurcation 
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(a) (b) 

Figure 4. Radially symmetric 2D system with rjkl = 0.2. In (a) is a bifurcation diagram for the system 
where dashed lines represent unstable branches and solid lines represent stable branches; numerical results are 
in black, and analytic results are in red (with circles). We find good agreement between the two for smaller 
e values, though the numerics display a large amplitude stable branch that the analytic solution does not. In (b) 
are plots of the numeric (black lines) and analytic (red with circles) solutions for Aam,p{t) with |e| = 0.01 and 
varying Q{0). The horizontal dashed line indicates the analytical unstable steady state for e < 0. The lower line 
corresponds to Q{0) = —0.01, e > 0, and the middle line corresponds to Q(0) = 0.3, e < 0; the agreement is 
good in these two cases. The upper line corresponds to Q(0) = 0.36, e < 0. This is above the unstable branch, 
so it grows to the large amplitude stable branch, which is not available from our analytical results. 


diagram for this geometry, which exhibits a transcritical bifurcation near the origin (black) 
that matches the theory (red with circles) well at small |e|. However, the numerics also 
display a large amplitude, stable bump solution that our theory does not predict. Unlike the 
subcritical ID case above, we do not extend to higher order approximations here. This large 
amplitude branch indicates that both the bump and ring steady state solutions are stable 
and available at positive e values, with the bump also being available in both a stable form 
and an unstable form over some range of negative e values. Figure 4(b) shows the evolution 
of both the numeric (black) and analytic (red with circles) Aamp(i) using |e| = 0.01 and 
three different values for <5(0). The first and lowest value is a ring with <5(0) = —0.01 (and 
positive e), which compares well with the analytic results. The second, intermediate value is a 
bump with Q(0) = 0.3 (and negative e), which is just slightly below the unstable branch, so we 
expect our analytic result above to work reasonably well in this case, and it does. However, 
the final value is a bump with Q(0) = 0.036 (and negative e), which is slightly above the 
unstable branch, so our analytic results above cannot be used. Numerically, though, we see 
that the solution grows until it reaches the stable, large amplitude branch. 

4.3. Fully 2D case. 

4.3.1. Squares. In this geometry, D is x G [0,vr/fc*], y G [0,vr/fc*]. The first case illus¬ 
trated is a supercritical system, in which = l/\/2 > 0.7 {rj = 0.01 and A;* = 10/2^/^^). The 
two plots in Figure 5 summarize the results here. Figure 5(a) shows the bifurcation diagram 
for this system at very small e values, and we find that there is good agreement at the smallest 
e values. Figure 5(b) plots the analytic and numerical solutions for ^amp(A) using e = 10“® 
and (5(0) = 10“^; there is very good agreement here, as e is very small in this case. 

The next case is a subcritical system, in which = 0.2 (rj = 0.01 and fe* = 2\/5), with 
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(a) (b) 

Figure 5. Fully 2D, square system with rjk^ = l/%/2. In (a) is a bifurcation diagram for the system where 
dashed lines represent unstable branches and solid lines represent stable branches; numerical results are in black, 
and analytic results are in red (with circles). We find good agreement for the smallest e values. In (b) are plots 
of the numeric (black) and analytic (red with circles) solutions for Aamp{t) with e = 10~® and Q(0) — 10“®; 
there is very good agreement here. 




(a) (b) 

Figure 6. Fully 2D, square system with gk^ = 0.2. In (a) is a bifurcation diagram for the system where 
dashed lines represent unstable branches and solid lines represent stable branches; numerical results are in 
black, and analytic results are in red (with circles). There is good agreement between the numeric and analytic 
solutions for smaller e values along the unstable branch, but the actual solution also displays a large amplitude, 
stable branch. In (b) are plots of the numeric (black) and analytic (red with circles) solutions for Aamp{t) with 
e = —0.001 and varying Q(0); the dashed line represents the analytic steady state value for this e. The lower 
line corresponds to Q(0) = 0.025, and there is good agreement here between the two. The upper line corresponds 
to Q(0) = 0.027, and the solution in this case rapidly grows to the large amplitude steady state. 


the results shown in Figure 6. Referring to Figure 6(a), the numerical solution (black) displays 
the small amplitude, unstable branch predicted by the theory above, and the numerics match 
the theory (red with circles) well at small e. However, there is also a stable, large amplitude 
branch in the bifurcation diagram that is not predicted by the theory, but which is similar to 
that seen in cases above. Figure 6(b) shows the evolution of Hamp(i) using e = —0.001 and 
two different values for Q{0). The first value is Q(0) = 0.025, which is just slightly below the 
unstable branch, so we expect our analytic results above (red with circles) to be close to the 
numerical results (black), and they are. However, the second <5(0) is 0.027, which is slightly 
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(a) (b) 

Figure 7. Fully 2D, hexagonal system with = 0.2. In (a) is a bifurcation diagram for the system 
where dashed lines represent unstable branches and solid lines represent stable branches; numerical results are 
in black, and analytic results are in red (with circles). We find good agreement between the two for smaller 
e values, though the numerics display both a large amplitude stable branch and secondary instabilities along the 
ring solution that the analytic solution does not. In (b) are plots of the numeric (black lines) and analytic (red 
with circles) solutions for Aamp(t) with |e| = 0.001 and varying Q{0). The horizontal dashed line indicates the 
analytical unstable steady state for e < 0. The lower line corresponds to Q(0) = —1.25 x 10““*, e > 0, with good 
agreement until t « 2500 (the ring solution is a saddle point in this geometry). The middle line corresponds 
to Q{0) = 6.125 X 10“®, e < 0, and the agreement is relatively good in this case. The upper line corresponds 
to Q(0) — 6.375 X 10“®, e < 0, which is above the unstable branch, so it grows to the large amplitude stable 
branch. 


above the unstable branch, so our analytic results cannot be used; we find here that the system 
rapidly approaches the large amplitude steady state seen in the bifurcation diagram. 

4.3.2. Hexagons. In this geometry, D is a: G [0,27r//i;*], y G [0,27r/\/3A;*]. As in the 
radially symmetric geometry, we have chosen a convention whereby bump solutions are shown 
with a positive Aamp value, while ring solutions are shown with a negative one. Simulations 
use rjkl = 0.2 {rj = 0.01 and /c* = 2\/5), and the results are shown in Figure 7. In Figure 7(a), 
we see our bifurcation diagram for this geometry, which exhibits a transcritical bifurcation 
near the origin (black) that matches the theory (red with circles) well at small |e|. However, 
the numerics display both a large amplitude stable bump solution and secondary bifurcations 
along the ring branch that our theory above cannot predict. As in the radially symmetric 
case, the large amplitude branch indicates that the bump solution is stable and available at all 
positive e values, in addition to some range of negative e values. The ring branch follows an 
interesting series of secondary bifurcations as e varies: the ring first breaks up into separate 
spots that grow in amplitude, the spots then begin moving toward each other to form a 
rectangle, and the spots of the rectangle eventually merge into the roll solution (the hnal 
continuous portion displayed for this branch is in fact just the roll solution). Of course, all 
of these solutions are unstable and, therefore, unlikely to be observed as steady states in any 
dynamical simulation. 

Figure 7(b) shows the evolution of both the numeric (black) and analytic (red with circles) 
Aamp(i) using |e| = 10“^ and three different values for (5(0)- The first and lowest value is 
a ring with <5(0) = —1.25 x (and positive e), which compares well with the analytic 
results until t ~ 2500. Due to the saddle point nature of the ring, the numerics eventually 
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begin to diverge from the analytic solution and end up at the large amplitude, stable state 
(albeit shifted from the standard bump solution due to the initial conditions). The second, 
intermediate value is a bump with <5(0) = 6.125 x 10“^ (and negative e), which is just slightly 
below the unstable branch, so we expect our analytic result above to work reasonably well 
in this case, and it does. However, the final value is a bump with Q(0) = 6.375 x 10“^ (and 
negative e), which is slightly above the unstable branch, so our analytic results above cannot 
be used. Numerically, though, we see that the solution grows until it reaches the stable, large 
amplitude branch. 

5. Hotspot suppression. Now that we know our system may exhibit two qualitatively 
different types of crime hotspots (supercritical and subcritical) it is natural to question what 
differences may exist, if any, between the behavior of these two classes of pattern with re¬ 
gards to hotspot suppression. As mentioned in the introduction, “hotspot policing” is a 
law-enforcement strategy whereby more police resources are focused on areas currently be¬ 
lieved to be within a hotspot in an effort to disrupt and destroy said hotspot. Field studies 
conducted to test the effectiveness of this strategy reveal that in some instances the hotspots 
seem to be destroyed, while in others they seem to simply be displaced. The 2D analyses 
we have performed above seem to offer an explanation as to why these two very different 
responses to suppression occur (refer to Figures 4 and 7). First, imagine a crime hotspot that 
exists within a linearly stable parameter regime (e < 0); the hotspot is therefore subcritical. 
If the police presence is enough to drive the attractiveness of the hotspot below the critical 
unstable branch of the bump solution in the bifurcation diagram, the system will tend to 
naturally drop down to the homogeneous state once suppression is relaxed, destroying the 
hotspot in question utterly. However, imagine now that the hotspot in question exists within 
the linearly unstable regime (e > 0) and is therefore supercritical. Any effort to suppress 
the bump solution will simply lead to the attractiveness being displaced to the surrounding 
area, i.e., a ring-like solution. Of course, in a fully 2D system, the ring solution will not be 
stable and will break up into separate spots, leading to a system that looks similar to the 
original one, but with the hotspots shifted to nearby regions. In this case, then, the hotspot 
policing will have simply led to a displacement of the hotspot to nearby areas rather than its 
destruction. 

The above hypothetical scenarios have been verified in computer simulations of the radially 
symmetric 2D continuous system and the fully 2D system for both the continuous and the 
discrete crime models. To do so, we choose a combination of parameters that are known to 
make the homogeneous steady state either linearly unstable or stable, whichever is desired. 
Then, we run simulations as described above in section 4 with initial conditions set to give a 
bump solution at the origin (considered to be the center of the field in the fully 2D case). We 
allow the simulation to run until a time ts when it seems to have reached a steady state, at 
which point we begin the suppression. This is accomplished by first defining an instantaneous 
damping field d(x) in the following way: 

(5.1) d(x) = ^ [1 - tanh [k (A(x, Q - Acutofr)]], 

where k sets the width of the transition region between total suppression and no suppression 
and Acutoff sets the attractiveness value above which suppression is desired. This damping 


Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 



478 


M. B. SHORT, A. L. BERTOZZI, AND P. J. BRANTINGHAM 




(b) 


Figure 8. Suppression in the radially symmetric 2D system with = 0.316, k = 3. The curves show 
A{r,t) as it evolves following the suppression that occurs at t = 0, and the horizontal dashed line represents 
Acutoff = 1. Shown in (a) is the case e = 0.4, and the suppression of the bump drives the system to the 
ring solution, which persists after the suppression is removed. Shown in (b) is the case e = —0.02, and the 
suppression of the bump drives the system to a temporary ring structure that decays to homogeneity once the 
suppression is removed. 


field is meant to represent police presence, which is concentrated almost exclusively in the 
areas of high attractiveness (hotspots). We assume that this presence has two effects. First, 
the damping held will reduce the crime rate in areas where there is a large police presence 
(d ~ 1). Second, the police presence will prevent burglars from beginning their search in these 
same areas. Mathematically, then, our PDF system is modihed to 


(5.2) 

1 

(M 

> 

II 

1 

(5.3) 

dp ^ 


2p, 

A 


-dpA + d(A- A°) 


Note that this damping held remains unchanged between any two successive values. In 
other words, the police may remain within an area for some time even after the crime there 
has been reduced. This is reasonable in the sense that in the real world, police do not have 
instantaneous information about what areas are most attractive and must instead rely on 
where events have occurred in the recent past when deciding where to allocate resources. 
Therefore, there is an inherent lag between the information possessed by the criminals and 
that possessed by the police. The typical timescale for this lag in the real world may be on 
the order of weeks to months [26], which is enough time for new hotspots to emerge [32, 5]. Of 
course, this damping method is only one of many possible choices, some of which are explored 
in [21], However, we suspect that for a large range of damping models, especially those with 
an appreciable temporal lag between criminal events and decisions on where police resources 
should be allocated, the basic outcomes described here will remain; subcritical spots may be 
destroyed, but supercritical spots will simply move. 

Results for the radially symmetric case are shown in Figure 8, and the hypothetical sce¬ 
narios play out as anticipated. In the supercritical case, suppression of the bump drives the 
system to the ring solution, which, due to its stability, remains even after suppression is 
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Figure 9. Suppression results for the fully 2D system with rjkl = 0.2, e = —0.05, k = 5, and Acutoff = 5.72. 
The top row shows results from the PDE system, while the bottom row shows results from the discrete system with 
equivalent parameters. These colormaps display high A in red, A in green, and low A in blue to purple. Shown 
in (a) is the system configuration right before suppression is first implemented. Soon after implementation, the 
central hotspot has disappeared entirely, but no further spots have emerged (b). Eventually the suppression is 
lifted, and the system begins to adopt the homogeneous steady state (c). 


relaxed (Figure 8(a)). Suppression of the subcritical bump initially sends the system to a 
ring-like state as well, since the suppression by definition will cause the origin to have very 
low A values, leaving the outer edge as the only place for criminal activity to occur. However, 
once the suppression is removed, the ring’s instability causes it to decay to the homogeneous 
state, and the original hotspot is now destroyed (Figure 8(b)). 

Figure 9 illustrates the effects of hotspot suppression in a fully 2D, subcritical system 
with periodic boundary conditions. Before suppression (Figure 9(a)), we see that our initial 
condition has led to a stable hotspot in the center of the field in both the continuum and 
discrete cases, though the discrete case also displays some quasi-hotspots near the edges of 
the domain due to random fluctuations that push the system at least temporarily above 
the unstable branch. Once suppression is introduced (Figure 9(b)), the hotspot dies away 
rather quickly, leaving an area of very low A in the center where the police presence remains 
and a faint ring near the domain edges. Critically, though, we do not see the emergence 
of new hotspots. Finally, when the next ts is reached (Figure 9(c)), there is actually no 
suppression needed since no hotspots remain, and the “coldspot” in the center returns to the 
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Figure 10. Suppression results for the fully 2D system with rjk^ = 0.2, e = 0.05, k = 5, and Acutoff = 6.12. 
The top row shows results from the PDE system, while the bottom row shows results from the discrete system with 
equivalent parameters. These colormaps display high A in red, A in green, and low A in blue to purple. Before 
suppression is first implemented, the system displays a number of hotspots (a). Soon after the implementation 
of suppression, the original hotspots vanish, but the attractiveness of the neighboring regions correspondingly 
increases, leading to a transient, ring-like structure that surrounds the location of the original central hotspot (b). 
By the time the next suppression time ts has arrived, a new steady state featuring hotspots near the original 
ones has been achieved (c). 


homogeneous value soon after the police leave the area. As predicted, the suppression was 
effective in eradicating the hotspot in the subcritical case. 

Figure 10 illustrates the effects of hotspot suppression in a fully 2D, supercritical system 
with periodic boundary conditions. Before suppression (Figure 10(a)), we see not only that 
our initial condition has led to a hotspot in the center of the field but that a number of 
other hotspots have developed near the edge due to the linear instability of the system. Once 
suppression is introduced (Figure 10(b)), the original hotspots disappear quickly. However, we 
see in the continuum case especially that the eradication of these spots has simply pushed the 
system into a different nonhomogeneous configuration, with a temporary ring-like structure 
surrounding the area where the central hotspot was located. Finally, by the time the next tg 
has arrived (Figure 10(c)), the system has reached a new steady state that exhibits hotspots 
in areas near where the original spots were. So, the suppression was ineffective in eradicating 
the supercritical hotspots and merely led to their displacement. 
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6. Conclusions. Through a weakly nonlinear analysis of our coupled PDE system ((2.4) 
and (2.5)), we have shown that in both the ID and 2D cases, our system may exhibit stable 
hotspots in both the supercritical and subcritical regimes. The existence of the subcritical 
hotspots offers another mechanism for crime pattern formation, in addition to the linear 
instability discussed in our previous work. 

Importantly, these distinct hotspot mechanisms may help explain the varying measures of 
success that police agencies have when attempting to suppress hotspots. In the supercritical 
case, suppression of a hotspot seems to simply displace the spot to neighboring regions, as the 
bump solution gives way to the ring solution, which either will be a new stable state (in the 
radially symmetric case) or will then break up into separate bumps (in the fully 2D case); this 
is illustrated in Figures 8(a) and 10. In the subcritical case, on the other hand, the suppression 
of the hotspot below the unstable bump solution branch of the bifurcation diagram (Figures 
4 and 7) should destroy it completely, as the only other stable state available in this regime 
is the homogeneous one; this is illustrated in Figures 8(b) and 9. 

As a corollary to this argument, we point out that the existence of these large amplitude 
branches introduces the possibility of hysteresis into the system. That is, if the parameters 
of the system are slowly varying with time (as social or economic conditions vary, perhaps), 
what was once a peaceful city (e < 0, homogeneous state) may experience a sudden burst of 
crime once the stability threshold is passed (as e > 0), rather than the crime slowly increasing 
as the parameters move further into the unstable regime (as would happen if the system were 
supercritical). Once in this linearly unstable state, police attempts at suppression may only 
have the effect of displacing crime hotspots. Furthermore, even if the parameters then change 
in such a way that e begins to decrease, this high level of crime may persist until the relevant 
parameters are even lower than when the initial outbreak occurred; though once the stability 
threshold is passed (so that e < 0), police suppression should help in destroying hotspots. 
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